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Abstract 


The function y = g(x) = \og(W(e x )), where W() denotes the 
Lambert W function, is the solution to the equation y + e y = x. 
It appears in various problem situations, for instance the calculation 
of current-volt age curves for solar cells. A direct calculation of g(x) 
may be inaccurate because of arithmetic underflow or overflow. We 
present a simple algorithm for calculating g{x ) that is robust, in that 
it will work for almost all x values which are representable in the 
arithmetic of one’s chosen computer language. The algorithm does not 
assume that the chosen computer language implements the Lambert 
W function. 


1 Introduction 


The Lambert W function w = f(z) is the solution of we w = 2 , for complex 
w and 0 . It can be considered as the multi-branch inverse of the conformal 
map w -A w e w = z between the complex w-plane and the complex 2 -plane. 
When w and 2 are restricted to having real values, the graph of the Lambert 
W function is as shown in figure [l] For further background regarding the 
Lambert W function, see DEI' 


Some problem situations, for instance the modeling of current and voltage in 
diodes or solar cells, reduce to an implicit equation which can be solved ex¬ 
plicitly by means of the Lambert W function. As a simple example, consider 
the implicit equation 

ae bU + bU = V (1) 
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where a, b are positive real numbers, parameters of the model, and U, V are 
real variables]^] The corresponding explicit solution for U as a function of V 

is 

^ = \ l°s(“ W(ae v ) S J (2) 

Here log() denotes the principal branch of the natural logarithm function, 
and W () denotes the principal branch of the Lambert W function. 


A typical task might be, given values of the model parameters, to draw a 
graph of U as a function of V. For that task, it is computationally efficient 
to have, instead of the implicit equation ([Tj) , the explicit solution (|2]) for U 
in terms of V and the model parameters. 


Another task might be to estimate the model parameter values a, b which 
best fit experimental observations of U, V pairs. For that task, one wants to 
have an understanding of how varying the model parameters will affect the 
U-V curves. 


Still another task might be to determine the relationships among the model 
parameters which correspond to having an extremum of a function of U 
or V. This task also requires an understanding of how varying the model 
parameters will affect the U-V curves, but it will also be helpful if the formula 
utilized has partial derivatives with respect to all the parameters. 


The expression in equation (j2j) is analytic, so is well-behaved with respect to 
its argument and the model parameters. It can be repeatedly differentiated, 
its extrema lie at stationary points, and so on. Because one is working with 
real values for U and V and the positive real model parameters, the argument 
for the Lambert W function evaluation in equation ([2]) is positive, so the 
principal branch of the Lambert W function is being used. Moreover, one is 
using the principal branch of the natural logarithm function. Everything is 
single valued in this expression. 


However, there can be a numerical difficulty in performing computations with 
equation (]2]) . Making the substitutions (coordinate changes) y = bU + log(a) 
and x = V + log(a), the computations involve an evaluation of the function 

V = 9{x) = \og{w {e x )^j . (3) 

2 Equation (|T|) is a considerable simplification, in order to have an example in mind. 
The typical solar cell model has four or five parameters. See [4], page 14, for instance. 
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Figure 1: Graph of y = LambertW (x) for both x and y real. 

The graph shows all pairs of real values of x, y 
which satisfy y — W(x)‘, that is, satisfy x = y e y 
Necessarily x > — 1/e. 

If — 1/e < x < 0, then there are two possible y values 
corresponding to the two branches of Lambert W function. 
The upper branch, with y > — 1, is the principal branch. 

If x is real and positive, the only real value of y = W(x) 
is the positive value of y on the principal branch. 


One difficulty which can arise, depending upon the computer programming 
language used, is numerical underflow or overflow, related to the evaluation of 
an exponential of x where x (negative or positive) has a large magnitude. The 
value of x must therefore be restricted to the set of logarithms of floating point 
numbers whose exponents can be accurately represented in the arithmetic 
facility of the computer language. A second difficulty which can arise is 
that the best computer programming language for the rest of one’s problem 
solution may not have a built-in Lambert W function evaluator. 


The purpose of this note is to address those two difficulties. We will describe 
a simple procedure, which can be implemented in any programming language 
with floating point arithmetic, for the robust calculation of the function 
y = log (IT (exp (a;))). The procedure is valid for essentially any real value 
of x which is representable in the programming language. 


2 Description of the Function y = g{x) = log (IT (exp (a;))) 


We may consider the function y = g{x) = log (IT (exp (a:))) as a transforma¬ 
tion of the Lambert W function, with a change of representation or coordinate 
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space. For clarity, we will restrict to real arguments. We can think of W(u), 
the principal branch of the Lambert W function, as a mapping of the posi¬ 
tive real line to itself. The function u = h(x) = exp(a;) maps the whole real 
line to the positive real line, and its inverse /r _1 (u) = log(w), the principal 
branch of the natural logarithm, maps the positive real line to the whole real 
line. In this interpretation, the function y = g(x) = log(W(exp(x))) is the 
composition of functions g = h^ 1 oW oh, and it maps the whole real line to 
the whole real line. 


Suppose that y = log(IF(exp(a:))). Taking exponentials (ie, applying the 
function h to both sides of the equation) gives 

e y = W(e x ). 

That is, using the definition of the Lambert W function, 


or 


e y+eV = e x . 


Taking logarithms (ie, applying the function h 1 to both sides) gives 


y + e y = x. 


( 4 ) 


Equation Q is a simple equation structure, as simple as the Lambert W 
defining equation structure 

w e w = z. (5) 

In fact, equation (|4]) is just equation (|5| in another coordinate system. 


When we are evaluating y = g(x) = log (W (exp(x))) as the solution to equa¬ 
tion Q, we are just evaluating the Lambert W function. There is an im¬ 
portant difference, however: The evaluation of g(x) does not involve much 
risk of underflow or overflow in the numerical representation of the computer 
language. 


Since y 


or 


g(x) satisfies g{x) + e 9 ^ = x, the first derivative of g(x) satisfies 
g\x) + e 9(a 4 g'{x) = 1 


g\x) = 


1 + e 9 ^ 


> 0 
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Figure 2: Graph of y — g{x) = log(fF(exp(;c))) 

for moderate magnitudes of the argument. 


The second derivative of g{x) satisfies 


g"{x) 


e 9 ( x ) g'(x) 

(1 + e ^ G)) 2 


< 0 


Figure [2] shows the function y = g(x ) for moderate values of the argument x, 
that is, for —4 < x < 4. Figure [3] shows the same function for larger values 
of the argument x, that is, for —1000 < x < 1000. 


One can see from these graphs that the function g(x) behaves like x when 
x is much less than 0, and behaves like log(x) when x is much more than 0. 
For values of x around 0, there is a smooth blend between the two behaviors, 
with g( 1) = 0. The function g{x) is strictly monotonic increasing, as it has 
a positive first derivative. It curves downward, as it has a negative second 
derivative. 


One can further see, from figure [3j that when the range of the argument x 
is large, the graph of y — g(x) looks like it has a sharp corner at the origin. 
Actually, as figure [2] illustrates, the graph does not really have a sharp corner. 
Nonetheless, at a suitable distance (large scale), one has in the function g 
a useful smooth function for representing a function which has a step in its 
derivative. 
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Figure 3: Graph of y = g{x) = log(lF(exp(a:))) 

for large magnitudes of the argument. 


3 Algorithm for Calculating y = g(x) = log (IF (exp (a;))) 


In the terminology of H. Kuki ([3j page 23), the function g(x ) is contracting. 
That is, \Ag(x)\ < | Arc |, or \g'(x)\ < 1, for all x values in its argument 
domain. That means the task of finding an estimate for g(x) given x is 
relatively stable. A slight change in x (noise in the input) will produce only a 
slight change in g(x). The only challenges in developing a formula to estimate 
g(x), given x, are finding an appropriate algorithm, coding the sequence of 
calculations to avoid unwanted cancellation, and being reasonably efficient 
in the number of computations performed. 


A suitable algorithm can be an initial estimate, followed by some number of 
iterations of a refinement. Halley’s method is used to perform refinements 
because it has cubic convergence, and the derivatives involved can be calcu¬ 
lated efficiently. 


Given any fixed real number x, we wish to find a real number y such that 

h(y) = y + e y — x 
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is zero. The first and second derivatives of h(y) are needed for Halley’s 
method. They are 

h\y) = 1 + e y 
h"{y) = e v 

and hence are particularly easy to calculate. Once one has e y from the 
calculation of h(y), the derivatives are also at hand. It is also necessary, in 
order to use Halley’s method, that the first derivative is non-zero; that is the 
case for h'(y) = 1 + e v . 

As an initial estimate y 0 , we choose to use y 0 — x for x < —e, and y 0 = 
log(x) for x > e. For — e < x < e, we linearly interpolate between the two 
values —e and 1. This is an extraordinarily crude initial estimate, but it is 
sufficient, since Halley’s method is very robust and rapidly convergent in this 
application. 


The general iteration formula for Halley’s method is 

= ,_ 2 h(y n )h'(y n ) 

y n +i Vn 2 h'(y n ) 2 - h{y n )h"{y n ) 

In this particular case h(y) = y + e y — x and the iteration formula becomes 


Vn+l Vn 


2 (y n + e Vn - x)(l + e Vn ) 

2(1 + e yn ) 2 — (y n + e Vn — x)e yn 


( 6 ) 


The details of coding depend upon the computer language. It will be efficient 
to evaluate e Vn only once per iteration. All other computations are straight¬ 
forward arithmetic. When evaluating the denominator of the adjustment in 
the iteration equation ([b]), there is little risk of cancellation resulting from the 
subtraction, as the first term in the denominator is larger than the second 
term. 

In practice, just a few iterations suffice to give a good result. For arguments 
in —10 6 < x < 10 6 , four iterations of Halley’s method reduce the absolute 
error to less than 10~ 80 . The actual coding can use a convergence criterion, 
based upon the desired maximum error in the estimate of function value, to 
determine how many iterations to perform. Alternatively, if the precision is 
fixed by the computer language’s arithmetic representation or by the needs 
of the application situation, then one can determine how many iterations of 
the refinement will suffice, and perform only that number, omitting the final 
redundant iteration which verifies the convergence. This technique, due to 
H. Kuki as seen in his algorithms for computing square root (see [3], pages 
49-50 and 135-136), probably deserves a name. Perhaps it should be called 
“Kuki’s convergence non-test”. 
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